CEATTLE is short for Climate-Enhanced, Age-based model with Temperature-specific Trophic Linkages and Energetics, which is a multi-species age-structured assessment model developed for groundfish in the Bering Sea, USA by Holsman et al. (2015). To incorporate the impacts of climate, the model includes temperature-dependent von Bertalanffy weight-at-age functions (VBGF) and temperature-specific, bioenergetics-based predation interactions. Inputs of the model include U.S. National Marine Fisheries Service Alaska Fisheries Science Center (AFSC) survey and fishery data. Outputs include historical estimates of predation mortality, fishing mortality, biomass, recruitment, etc.

Rceattle is an R package designed to implement the CEATTLE model using Template Model Builder (TMB; Kristensen et al. (2015)), which can be installed using following https://github.com/kaskr/adcomp/wiki/Download. Rceattle is structured similar to the original manuscript in terms of modularization. Seperate functions (i.e. modules) estimate retrospective temperature- and size-specific predator rations, prey preference, and weight-at- age. These are then used as inputs to the CEATTLE model to evaluate how predation mortality, recruitment, and survival of three target species change under historical climate conditions and harvest rates.

Currently RCeattle can take survey index, survey catch-at-age, fishery catch, fishery catch-at-age, bioenergetics, life history, and diet based inputs and estimates time series of biomass, fishing mortality, and predation mortality using the CEATTLE model in TMB. To run Rceattle a user can either load a pre-build Bering Sea data set, build a data set from an excel spreadsheet, or build a data set from ADMB CEATTLE dat files by selecting a control file located in a data directory with .dat files used in ADMB and using the build_dat function. The model will then initialize assuming $0$ for all parameters unless a parameter list is provided and/or "ceattle.par" or "ceattle_par.std" is selected. In the later case, a .par or .std file provided from a previously estimated ADMB model can be used and specified with associated directory (i.e. "mydir/ceattle_par.std"). Documentation for all functions can be found using ?.

Installation

To install Rceattle you will need to do the following steps:

# Install devtools if you don't already have it
install.packages("devtools")
# Install TMB and rtools
# Instructions can be found here for non-pc: https://github.com/kaskr/adcomp/wiki/Download
install.packages('TMB', type = 'source')
# Try "TMB::runExample(all = TRUE)" to see if TMB works

# Install Rceattle
devtools::install_github("grantdadams/Rceattle", auth_token = "4925b42ac46f1e0aefd671e9dc0c1cf1b3157017", build_vignettes = TRUE)

Data

For example, to run the 2017 single-species or multi-species assessment for the Bering Sea, a data file must first be loaded:

library(Rceattle)
data(BS2017SS) # ?BS2017SS for more information on the single-species data
data(BS2017MS) # ?BS2017MS for more information on the multi-species data

The data can then be written to excel, altered and reread into R in a useable format for Rceattle. The first page of the excel document will have meta-data describing each item.

# Write data to excel
Rceattle::write_excel(data_list = BS2017SS, file = "BS2017SS.xlsx")

# Change the data how you want in excel
# Read the data back in
mydata <- Rceattle::read_excel( file = "BS2017SS.xlsx")

Estimation

Single-species mode

Then the model can be fit by setting msmMode = 0 using the Rceattle function:

ss_run <- Rceattle::fit_mod(data_list = mydata,
                            inits = NULL, # Initial parameters = 0
                            file = NULL, # Don't save
                            debug = 0, # Estimate
                            random_rec = FALSE, # No random recruitment
                            msmMode = 0, # Single species mode
                            silent = TRUE)
# Type ?fit_mod for more details

The you can plot the model results using using

Rceattle::plot_biomass(Rceattle =  ss_run)
Rceattle::plot_recruitment(Rceattle =  ss_run)

Diagnostics

Survey diagnostics can be plot following:

Rceattle::plot_index(Rceattle =  ss_run)
Rceattle::plot_srv_comp(Rceattle =  ss_run)

Fishery diagnostics can be plot following:

Rceattle::plot_catch(Rceattle =  ss_run)
Rceattle::plot_fsh_comp(Rceattle =  ss_run)

Multi-species mode

For the a multispecies model starting from the single species parameters estimated, the following can be specified. Here I use the data in the single-species run, but the multi-species data can be loaded using data("BS2017MS"). Note: the only difference from BS2017SS is the residual mortality is lower. Note, starting with the single-species parameter estimates is ideal, it may not converge otherwise.

ms_run <- Rceattle::fit_mod(data_list = mydata,
                            inits = ss_run$estimated_params, # Initial parameters from single species ests
                            file = NULL, # Don't save
                            debug = 0, # Estimate
                            niter = 10, # 10 iterations around population and predation dynamics
                            random_rec = FALSE, # No random recruitment
                            msmMode = 1, # MSVPA based
                            suitMode = 0, # empirical suitability
                            silent = TRUE)

We can plot both runs as well:

mod_list <- list(ss_run, ms_run)
mod_names <- c("SS", "MS")

# Plot biomass trajectory
plot_biomass(Rceattle = mod_list, model_names = mod_names)
plot_recruitment(Rceattle = mod_list, model_names = mod_names, add_ci = TRUE)

plot_selectivity(Rceattle = mod_list, model_names = mod_names)

plot_mort(Rceattle = mod_list, model_names = mod_names, age = 2)

Projection

PROJECTION 1: CHANGING F Rceattle automatically projects the population forward when estimating. To check to see what the F rates are for each fishery we can check:

mydata$fsh_control$proj_F
mydata$projyr # Year the population is projected forward

# We can then change the F - For example, mean historical F
mydata$fsh_control$proj_F <- c(0.2342936, 0.513, 0.0774777)

# Re-run, without estimating
ms_run_proj <- Rceattle::fit_mod(data_list = mydata,
                            inits = ms_run$estimated_params, # Initial parameters from single species ests
                            file = NULL, # Don't save
                            debug = TRUE, # Do not estimate. Not changing parameters right now
                            niter = 10, # 10 iterations around population and predation dynamics
                            random_rec = FALSE, # No random recruitment
                            msmMode = 1, # MSVPA based
                            suitMode = 0, # empirical suitability
                            silent = TRUE)

PROJECTION 2: CHANGING FUTURE RECRUITMENT Change recruitment deviations - NOTE: the workflow for this may change in the future

# Get years from projection
nyrs <- mydata$endyr - mydata$styr + 1
nyrs_proj <- mydata$projyr - mydata$styr + 1
yrs_proj <- (nyrs + 1):nyrs_proj

# Replace future rec_devs with numbers
ms_run$estimated_params$rec_dev[,yrs_proj] <- replace(
  ms_run$estimated_params$rec_dev[,yrs_proj],
  values = rnorm( length(ms_run$estimated_params$rec_dev[,yrs_proj]),
                  mean = 0,
                  sd = 0.707) # Assumed value from penalized likelihood
  )

ms_run_proj2 <- Rceattle::fit_mod(data_list = mydata,
                                 inits = ms_run$estimated_params, # Initial parameters from single species ests
                                 file = NULL, # Don't save
                                 debug = TRUE, # Do not estimate. Not changing parameters right now
                                 niter = 10, # 10 iterations around population and predation dynamics
                                 random_rec = FALSE, # No random recruitment
                                 msmMode = 1, # MSVPA based
                                 suitMode = 0, # empirical suitability
                                 silent = TRUE)

# plot
mod_list <- list(ss_run, ms_run, ms_run_proj, ms_run_proj2)
mod_names <- c("SS", "MS no F", "MS mean historical F", "MS mean historical F w random rec")

# Plot biomass trajectory
plot_biomass(Rceattle = mod_list, model_names = mod_names, incl_proj = TRUE)
plot_recruitment(Rceattle = mod_list, model_names = mod_names, add_ci = TRUE, incl_proj = TRUE)

Retrospective analysis

We can also do retrospective peels and calculate Mohn’s Rho on the CEATTLE. NOTE: this is using mean historical F for projections as we changed it above

# Run retrospective
ms_run_proj_retro <- retrospective(Rceattle = ms_run_proj, peels = 10)

# Look at Mohns rho
ms_run_proj_retro$mohns

# Plot retrospectives
plot_biomass(ms_run_proj_retro$Rceattle_list)

# Seet how forecast changes
plot_biomass(ms_run_proj_retro$Rceattle_list, incl_proj = TRUE, mohns = ms_run_proj_retro$mohns)

Simulation

Data can be simulated from the estimated quantities using sim_mod:

# Simulate the single species model
ss_sim <- sim_mod(ss_run)

ss_sim_run <- Rceattle::fit_mod(
  data_list = ss_sim,
  inits = NULL, # Initial parameters = 0
  file = NULL, # Don't save
  debug = 0, # Estimate
  random_rec = FALSE, # No random recruitment
  msmMode = 0, # Single species mode
  silent = TRUE)

# Simulate the multispecies model
ms_sim <- sim_mod(ms_run)

ms_sim_run <- Rceattle::fit_mod(
  data_list = ms_sim,
  inits = ss_run$estimated_params, # Initial parameters = 0
  file = NULL, # Don't save
  debug = 0, # Estimate
  random_rec = FALSE, # No random recruitment
  msmMode = 1, # Holsman MS mode
  silent = TRUE)

# We can plot the simulated runs
mod_list <- list(ss_sim_run, ms_sim_run)
mod_names <- c("SS sim", "MS sim")

plot_biomass(Rceattle = mod_list, model_names = mod_names)
plot_recruitment(Rceattle = mod_list, model_names = mod_names)

Model variants

For recruitment, the model assumes that recruitment $R$ of species $sp$ follows the following equation:

R_{s p, y r} = R 0_{s p} * e^{R d e v_{s p, y r}}

$R_{sp,yr} = R0_{sp} * e^{Rdev_{sp, yr}}$

R d e v_{s p, y r} \sim N (0, σ_{r, s p})

$Rdev_{sp, yr} \sim N(0, \sigma_{r,sp})$ Following Holsman et al. (2015),

σ_{r, s p}

$\sigma_{r,sp}$ can be fixed at sqrt(0.5) and the model estimated using penalized likelihood, by setting ranfom_rec = FALSE. Alternatively, Rdev can be treated as random effects and

σ_{r, s p}

$\sigma_{r,sp}$ estimated by setting random_rec = TRUE:

# For recruitment, the model can estimate recruitment deviates as random effects
ss_re <- Rceattle::fit_mod(
  data_list = mydata,
  inits = NULL, # Initial parameters = 0
  file = NULL, # Don't save
  debug = 0, # Estimate
  random_rec = TRUE, # Turn of recruitment deviations as random effects
  msmMode = 0, # Single species mode
  silent = TRUE)

Diet uncertainty can also be added into by assuming a log-normal or gamma suitability function and multinomial error structure by changing suitMode. Type ?fit_mod for more information.

ms_gamma <- Rceattle::fit_mod(
  data_list = mydata,
  inits = ss_run$estimated_params, # Initial parameters from single species ests
  file = NULL, # Don't save
  debug = 0, # Estimate
  niter = 10, # 10 iterations around population and predation dynamics
  random_rec = FALSE, # No random recruitment
  msmMode = 1, # MSVPA based
  suitMode = 1, # Have a gamma function with time-independent length ratio for suitability. Includes diet in likelihood as multinomial
  silent = TRUE)
# Can try different functions. Look at ?fit_mod suitmode

References

Holsman, K.K., Ianelli, J., Aydin, K., Punt, A.E., Moffitt, E.A., 2015. A comparison of fisheries biological reference points estimated from temperature-specific multi-species and single-species climate-enhanced stock assessment models. Deep-Sea Research Part II: Topical Studies in Oceanography 134, 360–378. https://doi.org/10.1016/j.dsr2.2015.08.001

Kristensen, K., Nielsen, A., Berg, C.W., Skaug, H., Bell, B., 2015. TMB: Automatic Differentiation and Laplace Approximation. arXiv 70, 1–21. https://doi.org/10.18637/jss.v070.i05

Rceattle

Grant Adams, Kirstin Holsman